package tech.waterism.model;

import cn.hutool.core.io.resource.ResourceUtil;
import cn.hutool.core.util.ArrayUtil;
import tech.waterism.modelbase.PredictBase;
import tech.waterism.modelbase.PredictResult;
import com.alibaba.fastjson.JSON;
import lombok.Getter;
import lombok.Setter;
import org.apache.log4j.Logger;
import tech.waterism.model.dto.ModelOutputSimple;
import tech.waterism.model.util.ModelOutputHelper;
import tech.waterism.util.JacksonUtil;

import java.util.*;

/**
 * 三水源滞后演算模型（LAG_3）- 河道汇流
 *
 * <p>该模型用于三水源（地表径流、壤中流、地下水）的河道汇流计算，
 * 采用马斯京根法和滞后演算法相结合的方式进行河道演进。</p>
 *
 * <p><strong>模型特点：</strong></p>
 * <ul>
 *   <li>三水源分别演算：地表径流(RS)、壤中流(RI)、地下水(RG)</li>
 *   <li>考虑河道滞时效应(LAG)</li>
 *   <li>马斯京根河道演进</li>
 *   <li>支持多河段串联计算(MP)</li>
 * </ul>
 *
 * @author 初京刚
 * @version 2.0
 * @since 2023-06-23
 */
public class LAG_3 extends PredictBase {

    private static final Logger logger = Logger.getLogger(LAG_3.class);

    // ==================== 模型参数 ====================

    /**
     * 马斯京根法参数
     */
    private double C0;  // 马斯京根系数0
    private double C1;  // 马斯京根系数1
    private double C2;  // 马斯京根系数2

    /**
     * 壤中流消退系数（无量纲）
     * 取值范围: 0 < CI < 1
     */
    private double CI;

    /**
     * 地下水消退系数（无量纲）
     * 取值范围: 0 < CG < 1
     */
    private double CG;

    /**
     * 马斯京根权重系数（无量纲）
     * 取值范围: 0 < X < 0.5
     */
    private double X;

    /**
     * 河道演算系数（h）
     * 表示河道滞时参数
     */
    private double KK;

    /**
     * 滞时（时段数）
     * 表示水流在河道中传播的滞后时间
     */
    private double LAG;

    /**
     * 河道演算消退系数（无量纲）
     * 取值范围: 0 < CS < 1
     */
    private double CS;

    /**
     * 河段数量（无量纲）
     * 用于多河段串联计算
     */
    private int MP;

    /**
     * 流域面积 (km²)
     */
    private double F;

    /**
     * 计算时段长度 (h)
     */
    private int T;

    // ==================== 输入数据 ====================

    /**
     * 降水时间序列
     */
    @Setter
    @Getter
    private Date[] TM;

    /**
     * 三水源产流量 (mm)
     * R[0]: 地表径流 RS
     * R[1]: 壤中流 RI
     * R[2]: 地下水 RG
     */
    private double[][] R;

    /**
     * 初始河道流量序列 (m³/s)
     */
    private double[] QSIG0;

    /**
     * 初始河段流量 (m³/s)
     */
    private double[] QXSIG;

    /**
     * 预报开始时间
     */
    private Date startTm;

    /**
     * 预报开始数组下标
     */
    private int startIndex;

    /**
     * 初始地表径流 (m³/s)
     */
    private double QSP;

    /**
     * 初始壤中流 (m³/s)
     */
    private double QIP;

    /**
     * 初始地下水流量 (m³/s)
     */
    private double QGP;

    /**
     * 河道流量序列 (m³/s)
     */
    @Getter
    private double[] QSIG;

    /**
     * 河段流量序列 (m³/s)
     */
    private double[] QXSIGSS;
    
    /** 最后生成的标准化输出对象（用于接口返回） */
    private ModelOutputSimple lastOutput;

    /**
     * 构造函数，从参数Map中初始化模型
     *
     * @param data 模型参数Map
     */
    public LAG_3(Map<String, Object> data) {
        super();
        logger.info("初始化三水源滞后演算模型...");

        try {
            // 解析基本参数
            parseParameters(data);

            // 解析时间序列数据
            parseTimeSeriesData(data);

            // 解析初始状态
            parseInitialState(data);

            // 验证参数合理性
            validateParameters();

            logger.info("三水源滞后演算模型初始化完成");
            logger.debug("流域面积: " + F + " km², 滞时: " + LAG + " 时段, 河段数: " + MP);
        } catch (Exception e) {
            logger.error("三水源滞后演算模型初始化失败", e);
            throw new RuntimeException("模型初始化失败", e);
        }
    }

    /**
     * 解析模型参数
     * 
     * <p>解析以下关键参数：</p>
     * <ul>
     *   <li>CG - 地下水消退系数（日值，0-1）</li>
     *   <li>CI - 壤中流消退系数（日值，0-1）</li>
     *   <li>X - 马斯京根权重系数（0-0.5）</li>
     *   <li>KK - 河道演算系数/河道滞时参数 (h)</li>
     *   <li>LAG - 滞后时间（时段数）</li>
     *   <li>CS - 河道演算消退系数（0-1）</li>
     *   <li>MP - 河段数量（整数）</li>
     *   <li>F - 流域面积 (km²)</li>
     *   <li>clen - 计算时段长度 (h)</li>
     * </ul>
     *
     * @param data 参数Map
     */
    private void parseParameters(Map<String, Object> data) {
        logger.info("解析模型参数...");

        CG = parseDouble(data, "CG", "地下水消退系数");
        CI = parseDouble(data, "CI", "壤中流消退系数");
        X = parseDouble(data, "X", "马斯京根权重系数");
        KK = parseDouble(data, "KK", "河道演算系数");
        LAG = parseDouble(data, "LAG", "滞时");
        CS = parseDouble(data, "CS", "河道演算消退系数");
        MP = (int) parseDouble(data, "MP", "河段数量");
        F = parseDouble(data, "F", "流域面积");
        T = (int) parseDouble(data, "clen", "计算时段长度");

        logger.info("参数解析完成");
    }

    /**
     * 解析时间序列数据
     * 
     * <p>支持两种格式：</p>
     * <ol>
     *   <li>新格式：flowSeries - 包含time, rain, runoff, rSim的对象数组</li>
     *   <li>旧格式：dt + rSim - 分别存储时间和三水源径流的数组</li>
     * </ol>
     * 
     * <p>三水源径流数据结构：</p>
     * <pre>
     * R[0][i] - RS：第i时段的地表径流 (mm)
     * R[1][i] - RI：第i时段的壤中流 (mm)
     * R[2][i] - RG：第i时段的地下水 (mm)
     * </pre>
     *
     * @param data 数据Map
     */
    private void parseTimeSeriesData(Map<String, Object> data) {
        logger.info("解析时间序列数据...");

        // 1. 解析时间序列
        if (data.containsKey("flowSeries")) {
            // 新格式：flowSeries (time + rain + runoff + sources)
            parseFlowSeriesFormat(data);
        } else {
            // 旧格式：dt, rSim
            parseOldFormat(data);
        }

        // 2. 解析开始时间
        parseStartTime(data);

        logger.info("时间序列数据解析完成，序列长度: " + TM.length);
    }

    /**
     * 解析新格式的flowSeries数据
     * 
     * <p>flowSeries数据结构：</p>
     * <pre>
     * {
     *   "flowSeries": [
     *     {
     *       "time": "2025-06-05 08:00:00",
     *       "rain": 10.5,
     *       "runoff": 5.2,
     *       "rSim": [2.1, 1.8, 1.3]  // [RS, RI, RG]
     *     },
     *     ...
     *   ]
     * }
     * </pre>
     * 
     * <p>提取内容：</p>
     * <ul>
     *   <li>time - 时间序列TM[]</li>
     *   <li>rSim - 三水源径流R[3][]，其中R[0]=RS, R[1]=RI, R[2]=RG</li>
     * </ul>
     *
     * @param data 数据Map，包含flowSeries字段
     */
    private void parseFlowSeriesFormat(Map<String, Object> data) {
        logger.info("解析新格式 flowSeries 数据...");

        List<Map<String, Object>> flowSeriesList = (List<Map<String, Object>>) data.get("flowSeries");
        int length = flowSeriesList.size();

        TM = new Date[length];
        R = new double[3][length];  // R[0]=RS, R[1]=RI, R[2]=RG

        for (int i = 0; i < length; i++) {
            Map<String, Object> point = flowSeriesList.get(i);

            // 解析时间
            if (point.containsKey("time")) {
                String timeStr = (String) point.get("time");
                try {
                    TM[i] = new java.text.SimpleDateFormat("yyyy-MM-dd HH:mm:ss").parse(timeStr);
                } catch (Exception e) {
                    logger.error("时间解析失败: " + timeStr, e);
                    TM[i] = new Date();
                }
            }

            // 解析三水源径流 (从 rSim 数组中获取)
            if (point.containsKey("rSim")) {
                Object rSimObj = point.get("rSim");
                if (rSimObj instanceof List) {
                    List<Object> rSimList = (List<Object>) rSimObj;
                    if (rSimList.size() >= 3) {
                        R[0][i] = parseNumber(rSimList.get(0));  // RS - 地表径流
                        R[1][i] = parseNumber(rSimList.get(1));  // RI - 壤中流
                        R[2][i] = parseNumber(rSimList.get(2));  // RG - 地下水
                    }
                }
            }
        }

        logger.info("flowSeries 解析完成: " + length + " 个时段");
    }

    /**
     * 解析旧格式数据（dt + rSim）
     * 
     * <p>旧格式数据结构：</p>
     * <pre>
     * {
     *   "dt": ["2025-06-05 08:00:00", ...],
     *   "rSim": [[RS1, RS2, ...], [RI1, RI2, ...], [RG1, RG2, ...]]
     * }
     * </pre>
     * 
     * <p>兼容性说明：</p>
     * <ul>
     *   <li>支持JSON字符串格式（需先解析）</li>
     *   <li>支持JSONArray对象格式（直接读取）</li>
     * </ul>
     *
     * @param data 数据Map，包含dt和rSim字段
     */
    private void parseOldFormat(Map<String, Object> data) {
        logger.info("解析旧格式数据（dt + rSim）...");

        // 解析时间序列
        TM = parseArray(data, "dt", Date[].class, "降水时间");

        // 解析rSim数组 (三水源产流)
        Object rSimObj = data.get("rSim");
        if (rSimObj instanceof com.alibaba.fastjson.JSONArray) {
            com.alibaba.fastjson.JSONArray rSimArray = (com.alibaba.fastjson.JSONArray) rSimObj;
            R = new double[rSimArray.size()][];
            for (int i = 0; i < rSimArray.size(); i++) {
                com.alibaba.fastjson.JSONArray innerArray = rSimArray.getJSONArray(i);
                R[i] = new double[innerArray.size()];
                for (int j = 0; j < innerArray.size(); j++) {
                    R[i][j] = innerArray.getDoubleValue(j);
                }
            }
        } else if (rSimObj instanceof String) {
            // 如果是字符串，先解析为 JSON
            R = JSON.parseObject((String) rSimObj, double[][].class);
        } else {
            R = JSON.parseObject(JSON.toJSONString(rSimObj), double[][].class);
        }

        logger.info("旧格式数据解析完成");
    }

    /**
     * 解析预报开始时间
     * 
     * <p>支持多种格式：</p>
     * <ul>
     *   <li>字符串："2025-06-05 08:00:00"</li>
     *   <li>JSON数组：["2025-06-05 08:00:00"]</li>
     * </ul>
     * 
     * <p>功能：</p>
     * <ul>
     *   <li>确定预报计算的起始时刻</li>
     *   <li>计算起始时刻在时间序列中的索引位置</li>
     *   <li>用于区分历史模拟和未来预报</li>
     * </ul>
     *
     * @param data 数据Map，包含start字段
     */
    private void parseStartTime(Map<String, Object> data) {
        Object startObj = data.get("start");
        String startString;

        if (startObj instanceof String) {
            startString = (String) startObj;
        } else if (startObj instanceof com.alibaba.fastjson.JSONArray) {
            com.alibaba.fastjson.JSONArray startArray = (com.alibaba.fastjson.JSONArray) startObj;
            startString = startArray.getString(0);
        } else {
            startString = JSON.parseObject(JSON.toJSONString(data.get("start")), String.class);
        }

        try {
            java.text.SimpleDateFormat sdf = new java.text.SimpleDateFormat("yyyy-MM-dd HH:mm:ss");
            startTm = sdf.parse(startString);
        } catch (Exception e) {
            logger.warn("开始时间解析失败，使用第一个时间: " + e.getMessage());
            startTm = TM[0];
        }

        startIndex = ArrayUtil.indexOf(TM, startTm);
        logger.debug("预报开始时间: " + startTm + ", 开始索引: " + startIndex);
    }

    /**
     * 解析模型初始状态
     * 
     * <p>初始状态包括：</p>
     * <ol>
     *   <li>QSIG0 - 初始河道流量序列 (m³/s)</li>
     *   <ul>
     *     <li>用于考虑滞后效应LAG</li>
     *     <li>长度应 ≥ LAG，否则自动填充</li>
     *   </ul>
     *   <li>QSP - 初始地表径流 (m³/s)</li>
     *   <li>QIP - 初始壤中流 (m³/s)</li>
     *   <li>QGP - 初始地下水流量 (m³/s)</li>
     *   <li>QXSIG - 初始河段流量数组 (m³/s)</li>
     *   <ul>
     *     <li>长度为MP+1（河段数+1）</li>
     *     <li>用于马斯京根河道演算</li>
     *   </ul>
     * </ol>
     * 
     * <p>注意事项：</p>
     * <ul>
     *   <li>初始状态直接影响模拟结果的前期精度</li>
     *   <li>建议通过历史连续模拟获得合理的初始状态</li>
     * </ul>
     *
     * @param data 数据Map
     */
    private void parseInitialState(Map<String, Object> data) {
        logger.info("解析初始状态...");

        // 解析QSIG数组 (初始河道流量)
        QSIG0 = parseArrayFlexible(data, "QSIG", double[].class);

        // 解析初始三水源流量
        QSP = parseDouble(data, "QSP", "初始地表径流");
        QIP = parseDouble(data, "QIP", "初始壤中流");
        QGP = parseDouble(data, "QGP", "初始地下水流量");

        // 解析QXSIG数组 (初始河段流量)
        double[] qxsigArray = parseArrayFlexible(data, "QXSIG", double[].class);

        // 扩展QXSIG数组至MP+1长度
        QXSIG = new double[MP + 1];
        QXSIGSS = new double[MP + 1];
        for (int i = 0; i < MP + 1; i++) {
            if (i > qxsigArray.length - 1) {
                QXSIG[i] = QXSIG[i - 1];
            } else {
                QXSIG[i] = qxsigArray[i];
            }
            QXSIGSS[i] = QXSIG[i];
        }

        logger.info("初始状态解析完成");
    }

    /**
     * 验证参数合理性
     * 
     * <p>检查以下参数范围：</p>
     * <ul>
     *   <li>CI ∈ (0, 1) - 壤中流消退系数</li>
     *   <li>CG ∈ (0, 1) - 地下水消退系数</li>
     *   <li>CS ∈ (0, 1) - 河道演算消退系数</li>
     *   <li>X ∈ (0, 0.5) - 马斯京根权重系数</li>
     *   <li>LAG ≥ 0 - 滞时非负</li>
     *   <li>F > 0 - 流域面积大于0</li>
     * </ul>
     * 
     * <p>参数超出合理范围时记录警告，但不中断计算</p>
     */
    private void validateParameters() {
        if (CI <= 0 || CI >= 1) {
            logger.warn("壤中流消退系数CI超出合理范围 (0,1)，当前值: " + CI);
        }

        if (CG <= 0 || CG >= 1) {
            logger.warn("地下水消退系数CG超出合理范围 (0,1)，当前值: " + CG);
        }

        if (CS <= 0 || CS >= 1) {
            logger.warn("河道演算消退系数CS超出合理范围 (0,1)，当前值: " + CS);
        }

        if (X <= 0 || X >= 0.5) {
            logger.warn("马斯京根权重系数X超出合理范围 (0,0.5)，当前值: " + X);
        }

        if (LAG < 0) {
            logger.warn("滞时LAG不能为负数，当前值: " + LAG);
        }

        if (F <= 0) {
            logger.warn("流域面积F必须大于0，当前值: " + F);
        }
    }

    /**
     * 解析参数并记录日志
     */
    private double parseDouble(Map<String, Object> data, String key, String desc) {
        Object value = data.get(key);
        if (value == null) {
            throw new IllegalArgumentException("缺少必需字段: " + key + " (" + desc + ")");
        }
        double result = parseNumber(value);
        logger.debug(desc + " (" + key + "): " + result);
        return result;
    }

    /**
     * 将Object解析为double
     */
    private double parseNumber(Object value) {
        if (value instanceof Number) {
            return ((Number) value).doubleValue();
        } else {
            return Double.parseDouble(value.toString());
        }
    }

    /**
     * 解析数组字段 - 兼容字符串和直接的JSON数组
     */
    private <T> T parseArray(Map<String, Object> data, String key, Class<T> clazz, String desc) {
        Object value = data.get(key);
        T result;
        if (value instanceof String) {
            // 旧格式：字符串形式的JSON数组
            result = JSON.parseObject((String) value, clazz);
        } else {
            // 新格式：已经是JSON对象
            result = JSON.parseObject(JSON.toJSONString(value), clazz);
        }
        logger.debug(desc + " (" + key + "): 已解析");
        return result;
    }

    /**
     * 解析数组字段 - 支持JSONArray格式
     */
    private <T> T parseArrayFlexible(Map<String, Object> data, String key, Class<T> clazz) {
        Object obj = data.get(key);
        if (obj instanceof com.alibaba.fastjson.JSONArray) {
            com.alibaba.fastjson.JSONArray jsonArray = (com.alibaba.fastjson.JSONArray) obj;
            if (clazz == double[].class) {
                double[] array = new double[jsonArray.size()];
                for (int i = 0; i < jsonArray.size(); i++) {
                    array[i] = jsonArray.getDoubleValue(i);
                }
                return (T) array;
            }
        } else if (obj instanceof String) {
            return JSON.parseObject((String) obj, clazz);
        }
        return JSON.parseObject(JSON.toJSONString(obj), clazz);
    }

    /**
     * 测试入口 - 使用转换后的标准JSON文件
     */
    public static void main(String[] args) {
        logger.info("开始三水源滞后演算模型测试...");
        logger.info("读取转换后的JSON文件: json/LAG-input.json");

        // 读取转换后的标准JSON文件
        String content = ResourceUtil.readUtf8Str("json/LAG-input.json");
        Map<String, Object> data = JSON.parseObject(content, Map.class);

        LAG_3 model = new LAG_3(data);
        PredictResult result = model.predict();

        logger.info("模型计算完成");
        logger.info("输出结果: " + JSON.toJSONString(result));
        System.out.println(JSON.toJSONString(result));
    }

    @Override
    public void init() {
        // 初始化已在构造函数中完成
    }

    @Override
    public PredictResult predict() {
        logger.info("开始三水源滞后演算预测计算...");

        Map<String, Object> simulationResult = LAG3model();

        // ===== 传统数组格式输出（向后兼容）=====
        PredictResult result = new PredictResult();
        double[] QSim = JacksonUtil.convert(simulationResult.get("QSim"), double[].class);
        double[] RS = JacksonUtil.convert(simulationResult.get("RS"), double[].class);
        double[] RI = JacksonUtil.convert(simulationResult.get("RI"), double[].class);
        double[] RG = JacksonUtil.convert(simulationResult.get("RG"), double[].class);
        
        result.setQSim(QSim);
        result.setRSim(new double[][]{RS, RI, RG});
        
        // 时间序列
        result.setDt(TM);
        
        // 计算总径流深度（三水源之和）
        if (R != null && R.length == 3 && R[0] != null) {
            double[] runoffTotal = new double[R[0].length];
            for (int i = 0; i < R[0].length; i++) {
                runoffTotal[i] = R[0][i] + R[1][i] + R[2][i];
            }
            result.setRunoffSim(runoffTotal);
        }

        // 保存状态信息
        List<String> QXSIG = toStringList(JacksonUtil.convert(simulationResult.get("QXSIGS"), double[][].class));
        List<String> QSIG = toStringList(JacksonUtil.convert(simulationResult.get("QSIGS"), double[][].class));
        List<String> QSP = toStringList(JacksonUtil.convert(simulationResult.get("QSS"), double[].class));
        List<String> QIP = toStringList(JacksonUtil.convert(simulationResult.get("QIS"), double[].class));
        List<String> QGP = toStringList(JacksonUtil.convert(simulationResult.get("QGS"), double[].class));

        Map<String, List<String>> status = new HashMap<>();
        status.put("QXSIG", QXSIG);  // 河段下断面各时段初始流量
        status.put("QSIG", QSIG);    // 各时段初始总基流
        status.put("QSP", QSP);      // 地表坡面汇流
        status.put("QIP", QIP);      // 壤中流
        status.put("QGP", QGP);      // 地下水流
        result.setNewStatus(status);
        
        // ===== 新增：扁平化输出并保存到文件 =====
        try {
            ModelOutputHelper helper = new ModelOutputHelper("LAG_3", "1.0");
            
            // 设置时间序列和预热期（根据输入的 start 时间判断）
            helper.setTimeStamps(TM)
                  .setWarmUpPeriods(startIndex > 0 ? startIndex : 0);
            
            // 添加三水源流量（扁平化数组）
            helper.addDataArray("QS", RS)
                  .addDataArray("QI", RI)
                  .addDataArray("QG", RG)
                  .addDataArray("QSim", QSim);
            
            // 添加三水源径流深
            if (R != null && R.length == 3) {
                double[] runoffTotal = new double[R[0].length];
                for (int i = 0; i < R[0].length; i++) {
                    runoffTotal[i] = R[0][i] + R[1][i] + R[2][i];
                }
                helper.addDataArray("R", runoffTotal)
                      .addDataArray("RS", R[0])
                      .addDataArray("RI", R[1])
                      .addDataArray("RG", R[2]);
            }
            
            // 获取状态变量时间序列（包含初始状态）
            double[] QSS = JacksonUtil.convert(simulationResult.get("QSS"), double[].class);
            double[] QIS = JacksonUtil.convert(simulationResult.get("QIS"), double[].class);
            double[] QGS = JacksonUtil.convert(simulationResult.get("QGS"), double[].class);
            
            // 添加状态变量时间序列
            helper.addStateArray("QSP", QSS)
                  .addStateArray("QIP", QIS)
                  .addStateArray("QGP", QGS);
            
            // 设置最终状态
            Map<String, Object> finalState = new HashMap<>();
            finalState.put("QSP", QSS[QSS.length - 1]);
            finalState.put("QIP", QIS[QIS.length - 1]);
            finalState.put("QGP", QGS[QGS.length - 1]);
            helper.setFinalState(finalState);
            
            // 设置模型参数
            Map<String, Double> params = new HashMap<>();
            params.put("CS_径流消退系数", CS);
            params.put("CI_壤中流消退系数", CI);
            params.put("CG_地下径流消退系数", CG);
            params.put("LAG_滞时", LAG);
            params.put("F_流域面积", F);
            helper.setModelParameters(params);
            
            // 设置元数据
            Map<String, Object> metadata = new HashMap<>();
            metadata.put("area_km2", F);
            metadata.put("timestep_hours", T);
            metadata.put("lag_periods", LAG);
            metadata.put("reaches", MP);
            helper.setMetadata(metadata);
            
            // 构建并保存到文件
            ModelOutputSimple output = helper.build();
            this.lastOutput = output;  // 保存输出对象供接口使用
            String outputPath = "output/LAG_3_output.json";
            ModelOutputHelper.saveToFile(output, outputPath);
            logger.info("✅ LAG_3模型输出已保存到: " + outputPath);
        } catch (Exception e) {
            logger.error("保存模型输出时出错: " + e.getMessage(), e);
        }

        logger.info("三水源滞后演算预测计算完成，输出序列长度: " + result.getQSim().length);
        return result;
    }
    
    /**
     * 获取最后生成的标准化输出对象
     * 
     * <p>该方法用于接口返回，提供与output文件一致的格式</p>
     * 
     * @return ModelOutputSimple对象，包含完整的模型输出信息
     */
    public ModelOutputSimple getLastOutput() {
        return lastOutput;
    }

    /**
     * 三水源滞后演算核心算法 - LAG_3模型主计算函数
     *
     * <p>物理过程：</p>
     * <pre>
     * 产流输入（RS, RI, RG）→ 消退演算 → 滞后演算 → 马斯京根河道演进 → 河道流量
     * </pre>
     * 
     * <p>算法流程：</p>
     * <ol>
     *   <li>初始化河道流量序列QSIG[]</li>
     *   <li>计算单位转换系数 CP = F/T/3.6（径流深→流量）</li>
     *   <li>换算时段消退系数 CG, CI（日值→时段值）：</li>
     *   <pre>
     *   CG_时段 = CG_日^(T/24)
     *   CI_时段 = CI_日^(T/24)
     *   </pre>
     *   <li>计算马斯京根系数 C0, C1, C2</li>
     *   <li>逐时段计算：</li>
     *   <ul>
     *     <li>三水源消退演算：</li>
     *     <pre>
     *     QGP = QGP × CG + RG × (1-CG) × CP  （地下水）
     *     QIP = QIP × CI + RI × (1-CI) × CP  （壤中流）
     *     QSP = RS × CP                       （地表径流）
     *     </pre>
     *     <li>河道滞后演算：</li>
     *     <pre>
     *     QSIG1 = QSIG1 × CS + (QGP+QIP+QSP) × (1-CS)
     *     </pre>
     *     <li>马斯京根多河段串联演进</li>
     *     <li>陡降系数处理（针对2012年8月后的特殊情况）</li>
     *   </ul>
     *   <li>构建状态序列QSIGS[][]</li>
     * </ol>
     * 
     * <p>关键概念：</p>
     * <ul>
     *   <li><strong>消退演算</strong>：模拟产流量逐渐衰减的过程</li>
     *   <ul>
     *     <li>地下水（QG）：最慢，基流特征明显</li>
     *     <li>壤中流（QI）：中等速度</li>
     *     <li>地表径流（QS）：最快，直接响应</li>
     *   </ul>
     *   <li><strong>滞后演算</strong>：模拟河道水流传播的延迟效应</li>
     *   <ul>
     *     <li>LAG - 滞后时间（时段数）</li>
     *     <li>CS - 河道演算消退系数</li>
     *   </ul>
     *   <li><strong>马斯京根法</strong>：模拟河道水流演进</li>
     *   <ul>
     *     <li>考虑河道蓄水和流量传播</li>
     *     <li>支持多河段串联计算</li>
     *   </ul>
     * </ul>
     * 
     * <p>输出变量：</p>
     * <ul>
     *   <li>QS[] - 地表径流流量序列 (m³/s)</li>
     *   <li>QI[] - 壤中流流量序列 (m³/s)</li>
     *   <li>QG[] - 地下水流量序列 (m³/s)</li>
     *   <li>QSim[] - 总流量序列 (m³/s)</li>
     *   <li>QSS[], QIS[], QGS[] - 状态时间序列（包含初始状态）</li>
     *   <li>QSIGS[][] - 河道流量状态矩阵</li>
     *   <li>QXSIGS[][] - 河段流量状态矩阵</li>
     * </ul>
     *
     * @return 模拟结果Map，包含以下键值：
     *         <ul>
     *           <li>RS, RI, RG - 输入的三水源径流深 (mm)</li>
     *           <li>QS, QI, QG - 三水源流量 (m³/s)</li>
     *           <li>QSim - 总流量 (m³/s)</li>
     *           <li>QSS, QIS, QGS - 三水源状态序列 (m³/s)</li>
     *           <li>QSIGS - 河道流量状态矩阵</li>
     *           <li>QXSIGS - 河段流量状态矩阵</li>
     *         </ul>
     */
    public Map<String, Object> LAG3model() {
        logger.info("执行三水源滞后演算核心算法...");

        double[] RS = R[0];  // 地表径流
        double[] RI = R[1];  // 壤中流
        double[] RG = R[2];  // 地下水
        int Pnum = RS.length;

        // 初始化河道流量序列
        QSIG = new double[Pnum + (int) LAG];
        initializeQSIG();

        // 计算单位转换系数和演算参数
        double CP = F / T / 3.6;  // 径流深度(mm) -> 流量(m³/s) 转换系数
        CG = Math.pow(CG, T / 24.0);  // 时段消退系数
        CI = Math.pow(CI, T / 24.0);

        // 计算马斯京根系数
        calculateMuskingumCoefficients();

        // 初始化状态变量
        double[][] QSIGS = new double[Pnum + 1][QSIG0.length];
        double[][] QXSIGS = new double[Pnum + 1][MP + 1];
        double[] QG = new double[Pnum];
        double[] QS = new double[Pnum];
        double[] QI = new double[Pnum];
        double[] QGS = new double[Pnum + 1];
        double[] QSS = new double[Pnum + 1];
        double[] QIS = new double[Pnum + 1];

        // 设置初始值
        QSIGS[0] = Arrays.copyOf(QSIG0, QSIG0.length);
        QGS[0] = QGP;
        QSS[0] = QSP;
        QIS[0] = QIP;
        QXSIGS[0] = Arrays.copyOf(QXSIG, QXSIG.length);

        double QSIG1 = (LAG <= 1) ? QSIG[0] : QSIG[(int) (LAG - 1)];

        // 获取开始月份和年份（用于陡降处理）
        Calendar time0 = Calendar.getInstance();
        time0.setTime(TM[startIndex]);
        Integer Mth0 = time0.get(Calendar.MONTH) + 1;
        Integer Yth0 = time0.get(Calendar.YEAR);

        // 逐时段计算
        for (int J = 0; J < Pnum; J++) {
            // 三水源消退演算
            QGP = QGP * CG + RG[J] * (1.0 - CG) * CP;  // 地下水
            QIP = QIP * CI + RI[J] * (1.0 - CI) * CP;  // 壤中流
            QSP = RS[J] * CP;                           // 地表径流

            QG[J] = QGP;
            QS[J] = QSP;
            QI[J] = QIP;

            QGS[J + 1] = QGP;
            QSS[J + 1] = QSP;
            QIS[J + 1] = QIP;

            // 河道滞后演算
            QSIG1 = QSIG1 * CS + (QGP + QIP + QSP) * (1.0 - CS);
            double QTSIG = QSIG1;

            // 马斯京根河道演进
            QSIG[(int) (J + LAG)] = LCHCO(MP, QTSIG, QXSIGSS);
            QXSIGS[J + 1] = Arrays.copyOf(QXSIGSS, QXSIGSS.length);

            // 陡降处理（针对2012年8月后的特殊情况）
            Calendar time = Calendar.getInstance();
            time.setTime(TM[J]);
            Integer Mth = time.get(Calendar.MONTH) + 1;
            Integer Yth = time.get(Calendar.YEAR);

            if (Yth > 2012 && Mth > 8 && Yth0 > 2012 && Mth0 > 8) {
                QSIG[J] = QSIG[J] * 0.9;
                if (J % 100 == 0) {
                    logger.debug("应用陡降系数 0.9 于时段 " + J);
                }
            }
        }

        // 构建QSIGS状态序列
        for (int i = 1; i < QSIG.length - LAG; i++) {
            for (int j = 0; j < QSIG0.length; j++) {
                if (i < QSIG0.length) {
                    QSIGS[i][j] = (i + j < QSIG0.length) ? QSIG0[j + i] : QSIG[i + j];
                } else {
                    QSIGS[i][j] = (i + j < QSIG.length) ? QSIG[i + j] : QSIG[QSIG.length - 1];
                }
            }
        }

        // 封装结果
        Map<String, Object> simulationResult = new HashMap<>();
        simulationResult.put("RS", R[0]);
        simulationResult.put("RI", R[1]);
        simulationResult.put("RG", R[2]);
        simulationResult.put("QS", QS);
        simulationResult.put("QI", QI);
        simulationResult.put("QG", QG);
        simulationResult.put("QSim", QSIG);
        simulationResult.put("QGS", QGS);
        simulationResult.put("QSS", QSS);
        simulationResult.put("QIS", QIS);
        simulationResult.put("QSIGS", QSIGS);
        simulationResult.put("QXSIGS", QXSIGS);

        logger.info("三水源滞后演算核心算法执行完成");
        return simulationResult;
    }

    /**
     * 初始化QSIG河道流量数组
     * 
     * <p>QSIG数组用于存储河道出口流量时间序列，长度为 Pnum + LAG：</p>
     * <ul>
     *   <li>前LAG个元素：由初始状态QSIG0[]填充</li>
     *   <li>后Pnum个元素：由模型计算得到</li>
     * </ul>
     * 
     * <p>填充策略：</p>
     * <ul>
     *   <li>如果 LAG ≤ 1：直接使用QSIG0[0]</li>
     *   <li>如果 QSIG0.length ≥ LAG：从QSIG0[]复制前LAG个值</li>
     *   <li>如果 QSIG0.length < LAG：</li>
     *   <ul>
     *     <li>先复制QSIG0[]中的所有值</li>
     *     <li>不足部分用最后一个值填充（常数延拓）</li>
     *   </ul>
     * </ul>
     * 
     * <p>物理意义：</p>
     * <p>考虑滞后效应LAG，需要知道过去LAG个时段的流量，
     * 这些历史流量通过QSIG0[]提供</p>
     */
    private void initializeQSIG() {
        if (LAG <= 1) {
            QSIG[0] = QSIG0[0];
        } else {
            for (int i = 0; i < LAG; i++) {
                if (QSIG0.length >= LAG) {
                    QSIG[i] = QSIG0[i];
                } else {
                    if (i < QSIG0.length) {
                        QSIG[i] = QSIG0[i];
                    } else {
                        QSIG[i] = QSIG[i - 1];
                    }
                }
            }
        }
    }

    /**
     * 计算马斯京根法系数 C0, C1, C2
     * 
     * <p>马斯京根法原理：</p>
     * <p>基于河段蓄水量与上下断面流量的线性关系：</p>
     * <pre>
     * S = K[X·I + (1-X)·Q]
     * </pre>
     * 其中：
     * <ul>
     *   <li>S - 河段蓄水量</li>
     *   <li>I - 入流流量</li>
     *   <li>Q - 出流流量</li>
     *   <li>K - 流量传播时间（KK参数）</li>
     *   <li>X - 权重系数（0-0.5）</li>
     * </ul>
     * 
     * <p>差分方程：</p>
     * <pre>
     * Q₂ = C0·I₂ + C1·I₁ + C2·Q₁
     * </pre>
     * 
     * <p>系数计算公式：</p>
     * <pre>
     * FKT = KK - KK×X + 0.5×T
     * C0 = (0.5×T - KK×X) / FKT
     * C1 = (KK×X + 0.5×T) / FKT
     * C2 = (KK - KK×X - 0.5×T) / FKT
     * </pre>
     * 
     * <p>约束条件：</p>
     * <ul>
     *   <li>C0 + C1 + C2 = 1（质量守恒）</li>
     *   <li>C0, C1, C2 ≥ 0（稳定性条件，要求 0 ≤ X ≤ 0.5 且 Δt ≤ 2KX）</li>
     * </ul>
     * 
     * <p>参数说明：</p>
     * <ul>
     *   <li>KK - 河道演算系数/马斯京根参数K (h)</li>
     *   <li>X - 马斯京根权重系数（通常取0.2-0.3）</li>
     *   <li>T - 计算时段长度 (h)</li>
     * </ul>
     */
    private void calculateMuskingumCoefficients() {
        double FKT = KK - KK * X + 0.5 * T;
        C0 = (0.5 * T - KK * X) / FKT;
        C1 = (KK * X + 0.5 * T) / FKT;
        C2 = (KK - KK * X - 0.5 * T) / FKT;

        logger.debug("马斯京根系数 - C0: " + C0 + ", C1: " + C1 + ", C2: " + C2);
    }

    /**
     * 马斯京根河道演进算法 - 多河段串联计算（LCHCO）
     *
     * <p>算法原理：</p>
     * <p>将河道分成MP段，每段应用马斯京根法，
     * 上一段的出流作为下一段的入流，依次演算：</p>
     * <pre>
     * 河段1：Q₁ = C0·I + C1·Q₁,prev + C2·Q₁,out
     * 河段2：Q₂ = C0·Q₁ + C1·Q₂,prev + C2·Q₂,out
     * ...
     * 河段MP：Q_MP = C0·Q_{MP-1} + C1·Q_{MP,prev} + C2·Q_{MP,out}
     * </pre>
     * 
     * <p>计算步骤：</p>
     * <ol>
     *   <li>IM = MP + 1（总河段数）</li>
     *   <li>特殊情况：IM = 1时，直接返回入流（无演算）</li>
     *   <li>一般情况：逐河段应用马斯京根公式：</li>
     *   <ul>
     *     <li>Q1 = 当前入流RQ</li>
     *     <li>Q2 = 前一时段该河段的入流 QX[J-1]</li>
     *     <li>Q3 = 前一时段该河段的出流 QX[J]</li>
     *     <li>计算：RQ_new = C0×Q1 + C1×Q2 + C2×Q3</li>
     *     <li>更新状态：QX[J-1] = RQ（保存当前入流）</li>
     *     <li>RQ = RQ_new（作为下一河段的入流）</li>
     *   </ul>
     *   <li>更新最后河段出流：QX[IM-1] = RQ</li>
     * </ol>
     * 
     * <p>物理意义：</p>
     * <ul>
     *   <li>模拟洪水波在河道中的传播过程</li>
     *   <li>洪峰削减（由于河道蓄水）</li>
     *   <li>洪峰滞后（由于传播时间）</li>
     * </ul>
     * 
     * <p>状态变量QX[]说明：</p>
     * <ul>
     *   <li>QX[0] - 第1河段的前一时段入流</li>
     *   <li>QX[1] - 第2河段的前一时段入流</li>
     *   <li>...</li>
     *   <li>QX[MP] - 最后河段的出流（即总出流）</li>
     * </ul>
     *
     * @param MP 河段数量（MP+1为总河段数）
     * @param RQ 当前时段入流流量 (m³/s)
     * @param QX 河段流量状态数组（长度MP+1），存储各河段前一时段的入流/出流
     * @return 河道出口流量 (m³/s)
     */
    private double LCHCO(double MP, double RQ, double[] QX) {
        double IM = MP + 1;
        if (IM == 1.0) {
            QX[(int) (IM - 1)] = RQ;
        } else {
            for (int J = 1; J < IM; J++) {
                double Q1 = RQ;
                double Q2 = QX[J - 1];
                double Q3 = QX[J];
                QX[J - 1] = RQ;
                RQ = C0 * Q1 + C1 * Q2 + C2 * Q3;
            }
            QX[(int) (IM - 1)] = RQ;
        }
        return RQ;
    }

    /**
     * 将double数组转为String列表
     */
    private List<String> toStringList(double[] array0) {
        String[] array = new String[array0.length];
        for (int i = 0; i < array.length; i++) {
            array[i] = String.valueOf(array0[i]);
        }
        return Arrays.asList(array);
    }

    /**
     * 将double二维数组转为String列表
     */
    private List<String> toStringList(double[][] array0) {
        String[] array = new String[array0.length];
        for (int i = 0; i < array.length; i++) {
            array[i] = Arrays.toString(array0[i]);
        }
        return Arrays.asList(array);
    }
}

